CHC analysis
Dan Bruzzese
6/17/2020
using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.
MANOVA for species level
I cannot do site contrasts because of our uneven samplign scheme. The MANOVA throws out samples.
fit_full_species_man$pca_summary ## Importance of first k=7 (out of 35) components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.0736 0.7112 0.47335 0.37674 0.33980 0.30966 0.2955
## Proportion of Variance 0.4066 0.1784 0.07903 0.05006 0.04073 0.03382 0.0308
## Cumulative Proportion 0.4066 0.5850 0.66402 0.71409 0.75481 0.78864 0.8194summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 4 0.03980 40.088 28 776.62 < 2.2e-16 ***
## sex 1 0.81535 6.956 7 215.00 1.795e-07 ***
## hostRace:sex 4 0.79391 1.833 28 776.62 0.005633 **
## Residuals 221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PCA for Species Level
LDA for Species level
## Reference
## Prediction Cingulata Cornivora Mendax pom Zepheria
## Cingulata 25 0 0 0 0
## Cornivora 1 5 0 0 0
## Mendax 0 0 60 5 1
## pom 0 0 1 116 0
## Zepheria 0 0 0 1 16## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 9.610390e-01 9.390180e-01 9.273307e-01 9.820323e-01 5.281385e-01
## AccuracyPValue McnemarPValue
## 1.508463e-49 NaN
Bayesian Model
Collapsed all 35 chc measurements into a value column with a corresponding chc column detailing what chc was measured. Fitted a model to evaluate the efect of host and sex on chc values. The model accounts for uneven site effects and individual effects.
All Bayesian models were created in Stan computational framework (http://mc-stan.org/) accessed with brms package (Bürkner, 2017). To improve convergence and guard against overfitting, we specified mildly informative conservative priors. We ran 4 chains with 2000 interations. We checked for convergence and found all models mixed well. (Pareto K estimates are all good!)
#1 value ~ hostRace * sex + (1|site) +(1|chc) + (1|ID)
Model Converged
mod = readRDS("fit_chc_model1.rds")
plot(mod)
model fit is good = posterior values match real values
brms::pp_check(mod, resp ="value") #great## Using 10 posterior samples for ppc type 'dens_overlay' by default.
Model fit is good = r2 == 88%
bayes_R2(mod, summary = TRUE)## Estimate Est.Error Q2.5 Q97.5
## R2 0.8790622 0.0009515364 0.8771306 0.8808946model summary
tab_model(mod)| Â | value | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | -0.25 | -0.49 – -0.00 |
| hostRace: Cornivora | -0.13 | -0.29 – 0.03 |
| hostRace: Mendax | 0.17 | 0.09 – 0.26 |
| hostRace: Zepheria | 0.33 | 0.23 – 0.43 |
| hostRace: pom | 0.11 | 0.04 – 0.17 |
| sex: Male | 0.03 | -0.06 – 0.11 |
| hostRaceCornivora.sexMale | -0.14 | -0.35 – 0.06 |
| hostRaceMendax.sexMale | -0.07 | -0.17 – 0.03 |
| hostRaceZepheria.sexMale | -0.01 | -0.14 – 0.12 |
| hostRacepom.sexMale | -0.14 | -0.23 – -0.05 |
| Random Effects | ||
| σ2 | 0.44 | |
| Ï„00 | 0.07 | |
| ICC | 0.86 | |
| N site | 8 | |
| N ID | 231 | |
| N chc | 35 | |
| Observations | 8085 | |
| Marginal R2 / Conditional R2 | 0.021 / 0.879 | |
two way interactions
plot_model(mod, type = "pred", terms = c("hostRace", "sex"))
fixed effects
plot_model(mod, bpe = "mean", bpe.style = "dot", prob.inner = 0.4, prob.outer = 0.8)
Manova for pomonella
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 1 0.82295 3.3194 7 108 0.003071 **
## site 1 0.85799 2.5537 7 108 0.017937 *
## sex 1 0.52503 13.9576 7 108 8.859e-13 ***
## hostRace:site 1 0.69396 6.8041 7 108 1.085e-06 ***
## hostRace:sex 1 0.95188 0.7799 7 108 0.605479
## site:sex 1 0.93727 1.0326 7 108 0.412735
## hostRace:site:sex 1 0.90238 1.6691 7 108 0.124190
## Residuals 114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Pomonella
## Reference
## Prediction Apple Haw
## Apple 25 12
## Haw 14 71## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.786885246 0.503288443 0.703538181 0.855813669 0.680327869
## AccuracyPValue McnemarPValue
## 0.006229331 0.844519267
LDA predicting host and sex for pomonella
## Reference
## Prediction Apple_Female Apple_Male Haw_Female Haw_Male
## Apple_Female 23 2 9 1
## Apple_Male 1 1 1 2
## Haw_Female 8 0 32 3
## Haw_Male 1 3 5 30## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 7.049180e-01 5.725131e-01 6.155825e-01 7.840219e-01 3.852459e-01
## AccuracyPValue McnemarPValue
## 8.606186e-13 9.110307e-01
MANOVA for Zepheria
# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.74708 0.67708 5 10 0.6507
## site 1 0.83749 0.38808 5 10 0.8461
## Residuals 14LDA for Zepheria
## Reference
## Prediction Zepheria_Female Zepheria_Male
## Zepheria_Female 5 3
## Zepheria_Male 4 5## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.5882353 0.1793103 0.3292472 0.8155630 0.5294118
## AccuracyPValue McnemarPValue
## 0.4062810 1.0000000
MANOVA for Mendax
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.57727 6.1024 6 50 7.541e-05 ***
## site 2 0.10716 17.1232 12 100 < 2.2e-16 ***
## sex:site 2 0.74705 1.3082 12 100 0.2258
## Residuals 55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Mendax
## Reference
## Prediction Mendax_Female Mendax_Male
## Mendax_Female 20 15
## Mendax_Male 13 13## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.54098361 0.07072905 0.40849889 0.66935590 0.54098361
## AccuracyPValue McnemarPValue
## 0.55241652 0.85010674
MANOVA for cingulata
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.78806 1.2774 4 19 0.3136059
## site 1 0.38209 7.6817 4 19 0.0007371 ***
## sex:site 1 0.70133 2.0229 4 19 0.1319094
## Residuals 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Cingulata
## Reference
## Prediction Cingulata_Female Cingulata_Male
## Cingulata_Female 5 7
## Cingulata_Male 7 7## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.46153846 -0.08333333 0.26587122 0.66629178 0.53846154
## AccuracyPValue McnemarPValue
## 0.83731354 1.00000000
MANOVA for Cornivora
There is only one site
# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.54542 0.83346 2 2 0.5454
## Residuals 3LDA for cornivora
Cant show a plot
## Reference
## Prediction Cornivora_Female Cornivora_Male
## Cornivora_Female 1 0
## Cornivora_Male 1 3## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.8000000 0.5454545 0.2835821 0.9949492 0.6000000
## AccuracyPValue McnemarPValue
## 0.3369600 1.0000000